# https://chat.deepseek.com/a/chat/s/3b7c02b0-f3e6-4862-8169-5465ed3ba34d
# https://chat.deepseek.com/a/chat/s/8120474d-15bd-42df-984b-45d579a4b0b7
import numpy as np
import random

def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    # 修正1: 变量名修正和初始化
    dataMatrix = np.mat(dataMatIn)
    labelMat = np.mat(classLabels).transpose()
    b = 0
    m, n = np.shape(dataMatrix)
    alphas = np.mat(np.zeros((m, 1)))
    iter_count = 0  # 避免使用关键字iter
    
    # 修正2: 添加主循环
    while (iter_count < maxIter):
        alphaPairsChanged = 0
        for i in range(m):
            # 1、求出目标的预测值
            fXi = float(np.multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[i, :].T) + b)
            # 2、计算预测值和真实值之间的误差
            Ei = fXi - float(labelMat[i])
            
            # 3、KKT条件（拉格朗日的约束条件）检查，如果不满足约束条件，放弃当前值
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or \
               ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
                # 4、随机选择第二个alpha
                j = selectJrand(i, m)
                # 5、计算第二个目标的预测值
                fXj = float(np.multiply(alphas, labelMat).T * (dataMatrix * dataMatrix[j, :].T) + b)
                # 6、计算第二个目标的预测值误差
                Ej = fXj - float(labelMat[j])
                
                # 7、保存旧的值对
                alphaIold = alphas[i].copy()
                alphaJold = alphas[j].copy()
                
                # 8、计算L和H，L和H是为了约束alphas[j]，如果不满足[L,H]区间，则需要continue，放弃当前的值对
                if (labelMat[i] != labelMat[j]):
                    L = max(0, alphas[j] - alphas[i])
                    H = min(C, C + alphas[j] - alphas[i])
                else:
                    L = max(0, alphas[j] + alphas[i] - C)
                    H = min(C, alphas[j] + alphas[i])
                # 8、如果L=H，则需要continue，放弃当前的值对
                if L == H: 
                    print("L==H")
                    continue
                
                # 9、计算eta，eta是alphas[j]变化的步长（eta = 2 * K(x_i, x_j) - K(x_i, x_i) - K(x_j, x_j)，其中K是核函数（这里使用线性核，即点积））
                eta = 2.0 * dataMatrix[i, :]*dataMatrix[j, :].T - \
                      dataMatrix[i, :]*dataMatrix[i, :].T - \
                      dataMatrix[j, :]*dataMatrix[j, :].T
                # 10、eta反映了目标函数在优化方向上的曲率。如果eta为负，说明目标函数在该方向上是凸的，存在极大值点；如果eta为正，则无法保证极大值的存在，跳过当前优化对(i, j)，选择下一对。
                if eta >= 0:
                    print("eta>=0")
                    continue
                
                # 11、 用计算出来的eta更新alpha[j]
                alphas[j] -= labelMat[j]*(Ei - Ej)/eta
                alphas[j] = clipAlpha(alphas[j], H, L)
                
                # 12、 检查修正后的值，如果步长太小，直接跳过当前优化对(i, j)，选择下一对
                if (abs(alphas[j] - alphaJold) < 0.00001):
                    print("j not moving enough")
                    continue
                
                # 13、 更新alpha[i]
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])
                
                # 14、 更新b
                # 选择逻辑
                #如果 α_i ∈ (0, C)
                #样本 i 是支持向量，用 b1（最可靠）
                #如果 α_j ∈ (0, C)
                #样本 j 是支持向量，用 b2
                #否则
                #两者都在边界上（α=0 或 α=C），取平均值（更鲁棒）

                b1 = b - Ei - labelMat[i]*(alphas[i]-alphaIold) * \
                    dataMatrix[i, :]*dataMatrix[i, :].T - \
                    labelMat[j]*(alphas[j]-alphaJold) * \
                    dataMatrix[i, :]*dataMatrix[j, :].T
                
                b2 = b - Ej - labelMat[i]*(alphas[i]-alphaIold) * \
                    dataMatrix[i, :]*dataMatrix[j, :].T - \
                    labelMat[j]*(alphas[j]-alphaJold) * \
                    dataMatrix[j, :]*dataMatrix[j, :].T
                
                if (0 < alphas[i]) and (C > alphas[i]): 
                    b = b1
                elif (0 < alphas[j]) and (C > alphas[j]): 
                    b = b2
                else: 
                    b = (b1 + b2)/2.0
                
                alphaPairsChanged += 1
                print(f"iter: {iter_count} i:{i}, pairs changed {alphaPairsChanged}")
        
        # 修正13: 迭代计数逻辑
        if (alphaPairsChanged == 0): 
            iter_count += 1
        else: 
            iter_count = 0
    
    return b, alphas

# 辅助函数
def selectJrand(i, m):
    j = i
    while (j == i):
        j = int(random.uniform(0, m))
    return j

def clipAlpha(aj, H, L):
    if aj > H:
        aj = H
    if aj < L:
        aj = L
    return aj

# 数据加载函数
def loadDataSet():
    # 使用之前生成的30个样本点
    class1 = np.array([
        [1.2, 3.5], [1.8, 2.8], [2.2, 3.8], [1.5, 4.2], [0.8, 3.0],
        [2.0, 4.5], [0.5, 2.5], [1.0, 3.8], [2.5, 3.0], [3.0, 4.0],
        [0.7, 4.0], [1.2, 2.2], [2.8, 3.5], [3.2, 2.8], [4.0, 3.2]
    ])
    class2 = np.array([
        [3.0, 1.2], [3.5, 0.8], [4.2, 1.5], [2.8, 0.5], [4.0, 0.8],
        [3.8, 2.0], [5.0, 1.2], [4.5, 2.5], [5.5, 1.8], [6.0, 2.0],
        [2.5, 1.0], [3.2, 1.8], [4.8, 0.5], [5.2, 0.8], [6.5, 1.0]
    ])
    
    dataMat = np.vstack([class1, class2]).tolist()
    labelMat = [1.0]*15 + [-1.0]*15
    return dataMat, labelMat

# 使用示例
if __name__ == '__main__':
    dataArr, labelArr = loadDataSet()
    b, alphas = smoSimple(dataArr, labelArr, 0.6, 0.001, 40)
    print(f"偏置项b: {b[0,0]:.4f}")
    print(f"支持向量: {alphas }")
    
    # 计算权重向量w
    # 在SVM中，权重向量 w 是支持向量的线性组合:
    dataMatrix = np.mat(dataArr)
    w = np.zeros((2, 1))
    for i in range(len(alphas)):
        if alphas[i] > 0:
            w += np.multiply(alphas[i]*labelArr[i], dataMatrix[i, :].T)
    
    print("\n=== 训练结果 ===")
    print(f"分隔平面方程: {w[0,0]:.4f}x + {w[1,0]:.4f}y + {b[0,0]:.4f} = 0")
    print(f"支持向量数量: {np.sum(alphas > 0)}")
    print(f"偏置项b: {b[0,0]:.4f}")